Stability and dynamics of vortex clusters in nonrotated Bose-Einstein condensates 
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We study stationary clusters of vortices and antivortices in dilute pancake-shaped Bose-Einstein 
condensates confined in nonrotating harmonic traps. Previous theoretical results on the stability 
properties of these topologically nontrivial excited states are seemingly contradicting. We clarify this 
situation by a systematic stability analysis. The energetic and dynamic stability of the clusters is 
determined from the corresponding elementary excitation spectra obtained by solving the Bogoliubov 
equations. Furthermore, we study the temporal evolution of the dynamically unstable clusters. 
The stability of the clusters and the characteristics of their destabilizing modes only depend on 
the effective strength of the interactions between particles and the trap anisotropy. For certain 
values of these parameters, there exist several dynamical instabilities, but we show that there are 
also regions in which some of the clusters are dynamically stable. Moreover, we observe that the 
dynamical instability of the clusters does not always imply their structural instability, and that for 
some dynamically unstable states annihilation of the vortices is followed by their regeneration, and 
revival of the cluster. 



PACS numbers: 03.75.Lm, 03.75.Kk, 03.65.Ge 

I. INTRODUCTION 

Quantized vortices are topological defects in systems 
with a long-range quantum phase coherence. They have 
been extensively investigated in different systems and 
branches of physics such as helium superfluids [1, 2], 
superconductors [3], cosmology [4, 5], and optics [6]. 
Atomic Bose-Einstein condensates (BECs) are extremely 
convenient systems to investigate the characteristics of 
quantized vortices due to their experimental versatility; 
many properties of these systems can be manipulated 
with lasers and external magnetic fields. Furthermore, 
real-space imaging of BECs can be carried out by optical 
in situ or time of flight measurements. 

In addition to extensive theoretical analysis, vortices 
and vortex clusters composed of several vortices with the 
same topological charge have been realized and inves- 
tigated experimentally in gaseous BECs [7-10]. These 
vortex clusters are local minima of energy for rotated 
condensates and hence stable states. However, recently 
there has been theoretical suggestions that other kinds 
of vortex clusters can be stationary and stable states for 
nonrotated condensates [11, 12]. These clusters typically 
consist of vortices and antivortices in specific configura- 
tions such that the various forces acting on the vortices 
exactly balance each other. Based on the seeming robust- 
ness under external perturbations of some of these topo- 
logical excited collective states, they could also be con- 
sidered as solitonic states — indeed, they seem to be sta- 
bilized only by strong enough nonlincarity of the conden- 
sate. The development of techniques for vortex creation 
in dilute BECs, and especially phase-imprinting meth- 
ods [13], may enable direct construction of such vortex 
clusters in the future. Both from the experimental and 



theoretical points of view, their stability properties is a 
central issue. 

These vortex structures have been investigated in the 
nonintcracting limit [11, 14] as well as in interacting sys- 
tems [11, 12, 15, 16]. In the previous studies, three dif- 
ferent stationary vortex cluster configurations have been 
found in interacting pancake-shaped condensates. The 
so-called vortex dipole and quadrupole states, shown in 
Figs. 1(a) and 1(c), were originally introduced by Craso- 
van et al. [11, 12]. The stationary vortex quadrupole 
state was found to exist in both interacting and nonintcr- 
acting condensates, whereas the stationary vortex dipole 
state exists only in condensates interacting strongly 
enough. Furthermore, the authors of Refs. [11, 12] stud- 
ied the stability against small external perturbations of 
these configurations by integrating the Gross-Pitaevskii 
equation in real time, and concluded them to be stable 
and robust in the nonlinear regime. 

Subsequently, based on the Bogoliubov quasiparticlc 
spectra of the cluster states, the stability properties of the 
stationary clusters were studied by Mottonen et al. [16]. 
It was found that at least for the parameter values stud- 
ied, both the vortex dipole and vortex quadrupole states 
are energetically and dynamically unstable. Effects of 
the energetic instability can be neglected at low enough 
temperatures due to vanishingly small dissipation, but 
the dynamical instabilities can prevent these states from 
being long lived. In addition, another stationary cluster, 
a vortex tripole, was introduced — see Fig. 1(b). 

In order to clarify the seeming contradiction between 
the results of Refs. [11, 12] and Ref. [16], we present in 
this paper a systematic study of the stability properties of 
the above-mentioned stationary vortex clusters as func- 
tions of total particle number and trap anisotropy. Wc 
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FIG. 1: Density profiles of the condensate for the stationary 
vortex dipole (a), tripole (b) and quadrupole (c) states. Vor- 
tices and antivortices are denoted by plus and minus signs 
at the vortex cores, respectively. The interaction strength 
parameter is g = 170. The unit of length is the harmonic 
oscillator length in the x-direction. 



find the stationary vortex cluster states by directly mini- 
mizing an error functional for the Gross-Pitacvskii equa- 
tion, and compute the corresponding instability modes, 
i.e., the quasiparticle excitations responsible for the dy- 
namical instability, from the Bogoliubov equations. Fi- 
nally, we investigate the nature of the instability modes 
by computing the time development of slightly perturbed 
cluster states. It is found that the characteristics of the 
instability modes only depend on the effective strength of 
the particle interactions and the trap anisotropy. Some 
dynamical instability modes do not destroy the cluster 
structure and the condensate is therefore referred to as 
structurally stable, whereas some lead to annihilation 
and subsequent revival of the vortices in the cluster. 



II. MEAN FIELD THEORY 

The dynamics of weakly interacting gaseous Bose- 
Einstein condensates in the zero-temperature limit is ac- 
curately described by the Gross-Pitaevskii (GP) equation 
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where the nonlinear operator H. = H[ip] acting on the 
condensate order parameter ip is given by 
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The interaction parameter g is determined by the s-wave 
scattering length a as g — A-Kh 2 a/m. Above, m is the 
atomic mass, and the external trapping potential is de- 
noted by Vt r (r). The order parameter is normalized ac- 
cording to J \ip\ 2 dr = N, where N is the total number 
of atoms in the condensate. The total energy of the con- 
densate can be calculated as 



J dr V'*( r ) 
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The stationary states of the condensate with an eigen- 
value /x are solutions to the GP equation of the form 
■0(r, t) — e — '^/^(r). Energetic and dynamic stability 
of a given stationary state can be inferred from the cor- 
responding excitation spectrum given by the Bogoliubov 
equations 
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where C = Ti.+g\4'\ 2 — ft, the quasiparticle amplitudes are 
denoted by u q (r) and v q (r) , and uj q is the eigcnfrcqucncy 
corresponding to the quasiparticle state with the index q. 

If the quasiparticle spectrum contains excitations with 
positive norm J dr (|it g | 2 — \v q \ 2 ) > but negative eigen- 
frequency oj q , the corresponding stationary state is ener- 
getically unstable, whereas eigenfrequencies with nonzero 
imaginary part indicate dynamical instability. The occu- 
pation of dynamical instability modes increases exponen- 
tially in time, and hence small perturbations of a dynam- 
ically unstable stationary state typically result in large 
modifications in its structure. 

In pancake-shaped traps, the system in question has 
a 50(2) rotation symmetry. However, the vortex clus- 
ter states shown in Fig. 1 are not rotationally invariant. 
Thus, according to the Goldstone theorem [17, 18], there 
should exist collective low-frequency modes which tend 
to restore the broken symmetry by rigidly rotating the 
cluster states. These modes are intrinsic to the vortex 
cluster states, since they are absent, e.g., for an axially 
symmetric single vortex state. 

Existence of a symmetry breaking induced long-lived 
collective mode is manifested as a 1/fc 2 divergence in the 
response function of a symmetry restoring variable F [17] . 
Let us consider operators of the form [19] 
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where b q and b q are the creation and annihilation opera- 
tors for a quasiparticle with index q and /(r) is a complex 
function. In the presence of dynamical instability modes, 
we may generalize the result presented Ref. [19] to obtain 
the response function as 
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with A q = | j dr /*(r) {ip* (r)u q (r) + ip(r)v q (r)^ |" and 

B q = | /dr/(r)(<0*(r)Ug(r) + ip(r)v q (r))\ 2 . The Bohr 
frequency u> q is given by u) q — (E q —E )/h, where E q = 
Ke{huj q } is the energy of the excitation q and Eq is the 
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ground state energy. From Eq. (6) we observe that a nec- 
essary condition for the response function \ F t° diverge 
as l/u>, i.e., as 1/fc 2 , is the existence of excitations with 
Re{tL> 9 } = 0. In Sec. IV, we show that there indeed exist 
a zero-energy instability mode for all stationary vortex 
clusters and we therefore identify them as the Goldstonc 
modes corresponding to broken rotational symmetry. 



III. COMPUTATIONAL METHODS 

In Ref. [16], the stationary vortex clusters considered 
were found to be energetically unstable. Thus, they 
cannot be found with the usual methods based on en- 
ergy minimization. We employ the method introduced in 
Ref. [16], in which stationary solutions of the GP equa- 
tion (1) are found by minimizing the error functional 
F[ip,n] = J dr\(H[ip] - n)ip(r)\ 2 . The chemical poten- 
tial which minimizes the error functional is given by the 
familiar expression /i = j dril<*(r)li.[ip]ip(r)/N. Using 
the functional gradient 
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one obtains the stationary solutions, e.g., by the method 
of steepest descent or by directly solving the equation 
5F/8ijj* = 0. 

For computational simplicity, we consider a pancake- 
shaped system in the harmonic potential 
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with io x ,u y -c oj z . Especially, we assume the z- 
confinemcnt to be tight enough, such that the low en- 
ergy solutions of the GP equation and the Bogoliubov 
equations can be taken to be of the form a(z) = 
l/V2^Ie- z2 /( 2a ') in the z-dircction. If the harmonic 
oscillator length a z = sj h/(mui z ) in the axial direc- 
tion is much larger than the scattering length a, the 
condensate can be described by the usual GP equa- 
tion with g = Aith 2 a/m [20]. Substituting the ansatz 
t/j(r, t) = i/j 2D (x, y, t)a(z) into Eq. (1), one obtains the ef- 
fectively two-dimensional (2D) Gross-Pitaevskii equation 
in the dimcnsionlcss form 



l dr 2 



-v 2 + (i 2 + av) + m + — 



*l>, (9) 



where ip(x,y,t) = a x t/) 2D (x,y,t) and a x = y/H/(mu> x ). 
The dimcnsionlcss quantities denoted by the tilde are 
obtained from the original dimensional ones by scaling 
the length by a x , time by 1/lo x , and energy by hco x . 
The dimcnsionlcss constant A = uj v /u> x characterizes trap 
anisotropy in the xy-plane. The normalization condition 
J dxdy\ip\ 2 = 1 implies the effective dimensionless 2D 
coupling constant to be g = A^/TrNa/a z . Equation (9) 



implies that the only relevant parameters in the problem 
are the coupling constant g and the trap anisotropy A. 
Thus it suffices to investigate the cluster states as func- 
tions of these parameters. 

Based on the minimization of the error functional, we 
have computed stationary cluster states for different val- 
ues of the interaction strength g and trap asymmetry A. 
To obtain a proper starting point for the minimization of 
the error functional F[tp, fj], we first minimize the energy 
of ansatz wave functions for which the condensate phase 
is fixed. The found stationary states were computed 
to the relative final accuracy F[tp s , n}/(N[i 2 ) < 1CT 15 , 
which confirms that these states were, indeed, very ac- 
curate stationary states. The Bogoliubov equations were 
then solved for these states, identifying possible excita- 
tion modes with non-real eigcnfrcqucncics. Furthermore, 
we studied the nature of the found instability modes 
by solving the condensate temporal evolution for initial, 
slightly perturbed states of the form 
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where ip s (r) is the condensate order parameter for a sta- 
tionary state, and u q (r), v q (r) are the quasiparticle am- 
plitudes of an instability mode. The constant n q de- 
(7) scribes the initial population of the instability mode — by 



using the normalization J dr (\u q \ 2 + \v q \ 2 ) = 1 for the 
instability modes, we typically chose n q ~ 10~ 2 -10 -1 
corresponding to less than one percent of the atoms in 
the unstable mode. We also point out that contrary to 
the real-frequency modes, the energy of the perturbed 
stationary state ip(r, 0) in Eq. (10) is up to second order 
in n q equal to the energy of the stationary state i^ a (r) if 
the perturbation corresponds to a dynamical instability 
mode. Furthermore, our numerical computations have 
showed that in the presence of several instability modes 
in the spectrum of the state, it is typically sufficient to 
consider only the excitation with the largest |Im{w 9 }|. 
This is due to the fact that quasiparticle interactions nec- 
essarily excite also the fastest growing imaginary mode 
which seems to always dominate the condensate dynam- 
ics. The temporal evolution of the perturbed state was 
solved numerically from the time-dependent GP equation 
using a combination of a split-operator method and an 
implicit Crank-Nicolson scheme. 



IV. RESULTS 

The vortex cluster states were analyzed by varying the 
interaction parameter in the range < g < 300 for the ro- 
tationally symmetric trap with A = 1, and by separately 
varying the anisotropy parameter A for fixed interaction 
parameter values. Specifically, we studied the existence 
of the stationary vortex clusters for g = 170 (which was 
used in the previous studies in Ref. [16]) and the ex- 
istence of certain instability modes for a wide range of 
the interaction parameter. Excitation spectra for all the 
stationary clusters considered contained negative energy 
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excitations with positive norm, implying these configu- 
rations to be energetically unstable. On the other hand, 
dynamical instability of the clusters turned out to be a 
more subtle issue. The imaginary parts of the dynamical 
instability modes of the clusters in rotationally symmet- 
ric traps are presented in Fig. 2 as functions of the in- 
teraction strength g. Inspecting the temporal evolution 
of the cluster states in the presence of excited instability 
modes, we observe that there are two qualitatively dif- 
ferent dynamical instability modes for each cluster; one 
which turns out to lead to decay of the cluster and one 
which tends to rotate the cluster rigidly. We identify the 
latter modes as the Goldstone modes discussed in Sec. II, 
since they have zero energy and they vanish if the trap 
becomes sufficiently anisotropic. 
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FIG. 2: Absolute values of the imaginary parts of the eigen- 
frequencies corresponding to instability modes as functions of 
the interaction strength g for stationary dipole, tripole and 
quadrupole clusters confined in rotationally symmetric trap, 
with A = 1. Curves corresponding to the vortex dipole are 
marked with v (decay mode) and a (rotational mode), the 
vortex tripole with o (decay mode) and ■& (rotational mode) , 
and the vortex quadrupole with □ (decay mode) and o (ro- 
tational mode). The dashed and the dotted line indicate the 
lower limits of the interaction strength for which stationary 
vortex dipole and tripole were found to exist, respectively. 



A. Vortex dipole 

For a rotationally symmetric trap with A = 1, the sta- 
tionary vortex dipole state was found to exist only for 
interaction strengths g > 42, i.e., certain amount of non- 
linearity is required for the existence of the stationary 
dipole configuration. This agrees with the observations 
presented in Ref. [12]. For the vortex dipole, there exists 
two dynamical instability modes: one for all interaction 
strength values g ~ 42-300 considered (marked with A in 
Fig. 2), and one that occurs only in the range g ~ 50-80 
(marked with v in Fig. 2). The latter mode dominates 



in the region g ~ 50-80 due to its larger |Im{u> g }|. Time- 
development of slightly perturbed vortex dipole states 
shows that the mode existing outside the region g ~ 50- 
80 corresponds very accurately to rigid rotation of the 
dipole configuration, as shown in Fig. 3. The structure 
of the dipole remains intact in this rotation, and hence 
the stationary dipole state can be considered structurally 
stable, although it is dynamically unstable in this region. 
Thus dynamical instability does not necessarily imply the 
cluster structure to be unstable. 
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FIG. 3: Temporal evolution of the slightly perturbed station- 
ary vortex dipole state for g = 160 and A = 1. The interaction 
strength g — 160 is in the region where only the rotational 
mode exists. Time is denoted by t and it is given in the units 
of l/u) x . The structure of the cluster remains very close to 
the original stationary configuration. 

The nature of the instability mode dominating vortex 
dipole dynamics in the region g ~ 50-80 is very differ- 
ent from the rotational mode. Temporal evolution of a 
slightly perturbed dipole state in this regime is shown 
in Fig. 4. The dominant mode renders the vortices of 
the dipole first to annihilate each other, but eventually 
the dipole configuration reappears from the vortex-free 
state. For g ~ 50, the density distribution of this re- 
vived dipole configuration is almost the same as in the 
initial state, but the vorticity is changed such that the 
vortex becomes an antivortex and vice versa. Similar be- 
havior has been observed in the numerical simulations of 
light propagation in a graded-index medium, where the 
vortex dipole nested in a light beam undergoes period- 
ical collapse and revival along the direction of the light 
propagation [21]. The collapse and revival of the dipole 
configuration continues periodically for g ~ 50, whereas 
for increasing interaction strength, nonlincarity gradually 
prevents the topological excitation from recombining and 
eventually the whole excitation disappears for g > 80. 

In asymmetric traps, the stationary vortex dipole state 
corresponding to interaction parameter g = 170 was 
found to exist in the range A ~ 0.9-1.5. Interestingly, 
we found that even the tiny amount of asymmetry corre- 
sponding to A = 1.005 is sufficient to eliminate the rota- 
tional instability mode for all values of g. This is consis- 
tent with the fact that in an asymmetric trap, the system 
no longer has the SO (2) symmetry and the Goldstone 
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FIG. 4: Collapse and revival of the slightly perturbed sta- 
tionary vortex dipole state for g — 60 and A = 1. It should 
be noted that apart from oscillations of the vortex locations, 
the vortex dipole state remains essentially intact for a long 
time before it starts to decay. The vortices have opposite 
topological charges after the revival of the dipole state. 



of this mode excites also the dominant instability mode 
leading to the decay of tripole cluster to dipole cluster. 
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FIG. 5: Temporal evolution of the perturbed stationary vor- 
tex tripole state for g — 160 and A = 1. The vortex tripole 
starts to decay almost immediately to the vortex dipole state 
due to the strong dynamical instability. 



mode should vanish. Thus by using a slightly asymmetric 
trap the stationary vortex dipole state can be made fully 
dynamically stable in ranges g ~ 42-50 and g > 80. Ex- 
perimentally, the transfer between the different stability 
and instability regimes can be accomplished by adjusting 
the total particle number and the trapping frequencies. 

B. Vortex tripole 

The vortex tripole state exists only in rather strongly 
interacting condensates: in symmetric traps with A = 1, 
the stationary vortex tripole state was found to exist only 
for g > 108. The stationary vortex tripole has always two 
dynamical instability modes, as shown in Fig. 2. The 
real part of the cigcnfrcqucncy vanishes for both of these 
modes, but the imaginary part of the dominating mode 
is roughly an order of magnitude larger than that of the 
other one. The decay of the tripole configuration cor- 
responding to the dominating instability mode is shown 
in Fig. 5. Under small perturbations, one of the outer- 
most vortices starts to drift out of the condensate and 
eventually it reaches the surface of the cloud and excites 
surface modes which can be observed in the last panel in 
Fig. 5. The remaining two vortices form a topologically 
neutral vortex dipole. Furthermore, the stationary vor- 
tex tripole exists as a stationary state only for g > 108 
which lies in the structurally stable region of the sta- 
tionary vortex dipole. In the numerical simulations, the 
vortex tripole takes a time of order 1 /oj x for the popula- 
tion Kg 10 _2 -10 _1 of the excitation, before it starts to 
notably decay. Thus the decay of the vortex tripole to the 
vortex dipole is slow enough to be experimentally observ- 
able. The nature of the other instability mode, to which 
we refer as the slow mode, is more difficult to find out 
since this mode never seemed to dominate the temporal 
evolution. Numerical simulations showed that the slow 
mode tends to rotate the condensate, but the evolution 



In anisotropic traps, the vortex tripole state turned 
out to be stationary in the same range A ~ 0.9-1.5 for 
g w 170 as the vortex dipole. Anisotropy of the trap cor- 
responding to A > 1.01 removes the slow instability mode 
from the spectrum, but the dominating mode always per- 
sists. This observation together with the characteristics 
of the initial decay dynamics suggests that the slow mode 
is analogous to the rotational mode of the vortex dipole 
and quadrupole. Since the dominant mode persists in an 
asymmetric trap, the vortex tripole state cannot be made 
dynamically stable by tuning the interaction strength or 
the trap anisotropy. 

C. Vortex quadrupole 

The vortex quadrupole was found to be stationary 
for all interaction strengths for A = 1. These results 
arc in agreement with the previously reported calcula- 
tions of Crasovan et al. [11], in which stationary vor- 
tex quadrupoles were found to exist both in interact- 
ing and noninteracting condensates. In this respect, 
the quadrupole configuration differs essentially from the 
dipole and tripole clusters. The dynamical instability 
modes of the quadrupole cluster resemble those of the 
vortex dipole: The stationary vortex quadrupole state in 
rotationally symmetric trap has two instability modes, 
one which exist for all values of the interaction parame- 
ter, and one existing only in the range g ~ 50-280, where 
it dominates the condensate dynamics, see Fig. 2. The 
mode which exists for all g has Re{w g } = 0, whereas the 
other mode has energy Kc{fku q } ~ 10~ 1 hui r . 

In the region g ~ 50-280, the dominant mode of 
the vortex quadrupole drives the vortices to merge to- 
gether and annihilate each other, but eventually the 
quadrupole configuration reappear as shown in Fig. 6. 
This mode also generates oscillations with increasing am- 
plitude such that the condensate stretches and shrinks 
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as the annihilation-revival cycles proceed. The ampli- 
tude of these oscillations increases gradually, and finally 
the vortices are driven out of the condensate. In the 
regions g < 60 and g ~ 280-300, the existing instability 
mode is analogous to the dipole cluster mode which corre- 
sponds to rigid rotation of the cluster, i.e., the Goldstonc 
mode. Thus the quadrupole cluster is structurally stable 
in this regime. Furthermore, one observes from Fig. 2 
that the vortex quadrupole tends to become dynamically 
stable in the limit of noninteracting condensate. It has, 
however, been shown that in the noninteracting case the 
persistent current such as the vortex quadrupole is al- 
ways structurally unstable against perturbations in the 
external trap parameters [22, 23]. Thus dynamical sta- 
bility does not necessarily imply stability against small 
perturbations in the trap parameters and vice versa, dy- 
namically unstable states can be structurally robust with 
respect to small perturbations. 
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FIG. 6: Collapse and revival of the vortex quadrupole for 
g = 260 and A = 1. In contrast to the vortex dipole state, the 
vortices have the same topological charge before and after the 
collapse and revival of the cluster. 



tion between vortices with same topological charge. 

Analogously to the vortex dipole case, computations 
with different interaction strength values g showed that 
for the stationary vortex quadrupole state, even the 
rather small amount of asymmetry A = 1.1 is sufficient to 
remove the rotational instability mode from the quasipar- 
ticle spectrum. Hence, also the vortex quadrupole state 
can be made dynamically stable in the regimes g < 60 
and g > 280. 

V. CONCLUSIONS 

We have systematically studied the existence, stabil- 
ity, and dynamics of stationary vortex clusters in dilute 
pancake-shaped BECs confined in nonrotating harmonic 
traps. In contrast to the previous investigations [11, 12], 
our approach utilizes the Bogoliubov equations to de- 
termine the quasiparticle spectrum which reveals unam- 
biguously not only the energetic and dynamic stability 
of a given state, but also yields explicitly the quasipar- 
ticle amplitudes of the instability modes. The nature of 
these modes was found out by computing the temporal 
evolution of slightly perturbed stationary cluster states. 

It was observed that the dipole, tripole and quadrupole 
clusters have various regimes of dynamical instability, but 
in some of these regimes the cluster configurations are 
structurally stable as they only tend to rotate rigidly. 
On the other hand, the vortex annihilation modes of the 
vortex dipole and quadrupole break the cluster structures 
altogether, but if energy dissipation is negligible, the clus- 
ters can reappear after a certain time. Furthermore, it 
was observed that even very small trap anisotropics suf- 
fice to remove the rotational instability modes, thus dy- 
namically stabilizing the dipole and quadrupole clusters 
for suitable values of the interaction parameter. 



As a function of the trap anisotropy parameter, the 
quadrupole cluster with g = 170 was found to be sta- 
tionary in the range A ~ 0.9-1.2. For increasing trap 
anisotropy parameter, the two vortices of the vortex 
quadrupole in the direction of the tight confinement move 
away from the center of the cloud, and eventually the 
state cannot be considered as a vortex quadrupole. This 
behavior is due to the increase of the buoyancy force in 
the direction of tight confinement and repulsive intcrac- 
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